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Flow of a thin viscous film down a flat inclined plane becomes unstable to long wave 
interfacial fluctuations when the Reynolds number based on the mean film thickness be¬ 
comes larger than a critical value (this value decreases as the angle of inclination with the 
horizontal increases, and in particular becomes zero when the plate is vertical). Control of 
these interfacial instabilities is relevant to a wide range of industrial applications includ¬ 
ing coating processes and heat or mass transfer systems. This study considers the effect 
of blowing and suction through the substrate in order to construct from first principles 
physically realistic models that can be used for detailed passive and active control studies 
of direct relevance to possible experiments. Two different long-wave, thin-film equations 
are derived to describe this system; these include the imposed blowing/suction as well 
as inertia, surface tension, gravity and viscosity. The case of spatially periodic blowing 
and suction is considered in detail and the bifurcation structure of forced steady states is 
explored numerically to predict that steady states cease to exist for sufficiently large suc¬ 
tion speeds since the film locally thins to zero thickness giving way to dry patches on the 
substrate. The linear stability of the resulting nonuniform steady states is investigated 
for perturbations of arbitrary wavelengths, and any instabilities are followed into the 
fully nonlinear regime using time-dependent computations. The case of small amplitude 
blowing/suction is studied analytically both for steady states and their stability. Finally, 
the transition between travelling waves and non-uniform steady states is explored as the 
suction amplitude increases. 


1. Introduction 

The flow of a viscous liquid film down an inclined plane, under the action of gravity, 
inertia and surface tension, is a fundamental problem in fluid mechanics that has received 
considerable attention both theoretically and experimentally due to the richness of its 
dynamics and its wide technological applications, e.g. in coating processes and heat or 
mass transfer enhancement. For sufficiently thin films (with other parameters such as 


inclination angle and viscosity fixed), the uniform Nusselt solution (Nusselt 1916) is 
stable, but for thicker layers the flow is susceptible to interfacial instabilities in the form 
of two-dimensional travelling waves which propagate down the slope, followed by more 
complicated time-dependent and three-dimensional behaviour. The linear stability of a 


uniform film was first considered by Benjamin (1957) and Yih (1963), who used an Orr 


Sommerfeld analysis to show that instability first appears at wavelengths that are large 
compared to the undisturbed film thickness h s . Using the Nusselt velocity at the free 
surface we define a Reynolds number R = p 2 gh? s sin#/2/x 2 (see (2.41 also) where p is 


the fluid density, g is the gravitational acceleration, 9 is the angle of inclination to the 
horizontal (see Figure [l]), and p is the viscosity of the fluid; the flow becomes linearly 
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unstable to long waves when R > R c = (5/4) cot 9 and we can see that the critical 
Reynolds number tends to zero as the plate becomes vertical. This result also shows that 
for a given angle the flow can be made unstable by increasing the density and/or the film 
thickness, or decreasing the viscosity. 

In order to go beyond linear theory without recourse to direct numerical simulations 
of the Navier-Stokes equations, a hierarchy of long-wave reduced-dimension models have 
been developed to analyse in detail the stability and nonlinear development of the flow 


(see reviews by 

Craster & Matar 2009 

Kalliadasis et al. 2012 

. The simplest fully nonlin- 

ear long-wave model was developed by 

Benney 

(1966 

) who used an expansion in a small 


slenderness parameter to derive a single evolution equation for the interface height. The 
Benney equation is valid for Reynolds numbers near the critical value R c (in fact it cap¬ 
tures exactly the linear stability threshold), in the sense that far away from critical in 
the unstable regime solutions can become unbounded in finite time ( Pumir et a/.| [l983), 
a phenomenon that invalidates the long wave approximation and is not observed in nu¬ 
merical simulations of the Navier-Stokes equations (Oron & Gottlieb 2002; Scheid et al. 
[2004). The Benney equation forms a rational basis for the development of asymptotically 
correct weakly nonlinear models that lead to the Kuramoto-Sivashinsky (KS) equation 
iTseluiko & Papageorgiou 2006, and references therein, for example). The KS equa- 


(see 


tion displays very rich dynamical behaviour including spatiotemporal chaos. In other 
canonical asymptotic weakly nonlinear regimes one can derive the generalised (i.e. dis- 
persively modified) KS equation along with electric-field induced instabilities (Tseluiko 
|fc Papageorgiou||2006 2010). In order to overcome the near-critical restrictions of the 
Benney equation, |Shkadov (1969) developed a coupled fully nonlinear long-wave system 
for the free-surface height and the local mass flux. The Slikadov model avoids finite-time 
singularities, but underpredicts the critical Reynolds number. Using a weighted-residual 


method, Ruyer-Quil & Manneville (2000) recently developed a new two-degree of free¬ 


dom long-wave model, which has the correct stability threshold and is well behaved in 
the nonlinear regime. In this paper we consider two such reduced-dimensional models in 
the presence of blowing and suction. 

In applications it is useful to be able to control the film dynamics. For instance in 
coating applications a stable uniform film is required, whereas in heat or mass transfer 
applications efficiency is improved if the flow is nonuniform and attains increased sur¬ 
face area and recirculation regions. Such diverse requirements motivate the introduction 
of extrinsic modifications to the system in the interest of controlling the dynamics. An 
example of such a modification is the utilization of a heated substrate which affects the 


interfacial dynamics via a combination of Marangoni effects and evaporation (Kalliada- 


sis et al. 2012). Substrate heating thus introduces new modes of instability relating to 


convection, and lowers the critical Reynolds number. The substrate behaviour can also 
be altered by allowing chemical coatings, elastic deformations, or interactions with flow 
through porous media ( Thiele et aL|2009 Ogden et aL|20lf Samanta et a/.|2()TT 2013), 
which is often modelled by an effective slip condition. External fields can also be used 
to stabilise or destabilise the interface. Depending on the fluids used, an applied mag¬ 
netic field can stabilise the flow (see Hsieh|l965 Shen et a77||199 Tj [Renardy fe Sun 1994]), 
whereas an electric field applied normal to the interface can destabilize the flow and in 
fact drive it to chaotic spatiotemporal dynamics even below critical R < R c (Tseluiko & 
Papageorgiou||2006 2010). 


One way to modify the dynamics of falling film flows is to topographically structure the 
substrate. There have been several theoretical and experimental studies of film flows down 
wavy inclined planes (typically with sinusoidal and step topographies), aiming to explore 
how topography affects stability and stability criteria such as critical Reynolds numbers, 
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how substrate heterogeneity interacts with nonlinear coherent structures, and from a 
practical perspective how topography induces flows that can be useful in heat or mass 


transfer by creating regions of recirculating fluid (see for example Tseluiko et al. 2013 and 


numerous references therein). The problem is quite complex with several parameters and 
an overall conclusion of these studies is that topography can either decrease or increase 
the critical Reynolds number in different regimes. 

Inhomogeneous heating of the substrate can generate nonuniform film profiles even in 
the absence of inclination (Sapryki n et a(.||2007 ). Scheid et al. (2002) used an extension 
of the Benney model to analyse flow over a substrate with a sinusoidally varying heat 
distribution, where temperature is coupled to flow via Marangoni effects. They solved 
the equations to obtain travelling waves in the case of uniform heating, and steady non- 
uniform solutions in the case of non-uniform heating, and used initial value calculations 
to demonstrate that imposing heating is able to halt the progress of a travelling wave, 
leading to stable steady interface shapes. The combination of localized heating and to¬ 
pography has also been considered and it has been demonstrated that features that form 
in the isothermal case (e.g. ridge formation ahead of step-down topography) can be re¬ 
moved by suitable heating (see Blyth fe Bassom||2012 Ogden et a/.||20ll). The removal 


of such a ridge has also been shown to be possible by the imposition of vertical electric 
fields rather than heating, providing another physical mechanism for interface control 
(see Tselu iko et q/.| [2b 08a| 6). 

Suction and injection blowing through an otherwise rigid substrate has well known 
applications in stabilizing flows and changing global structures that can negatively affect 
performance, such as boundary layer separation for example. In the types of interfacial 
flows of interest here there has been much more limited exploration; |Momoniat et ai\ 
(2010) studied the effect of imposing either suction or injection on a spreading drop and 


found that injection enhances ridge formation, while suction leads to cavities forming on 
the free surface. As the total mass is not conserved, a steady state is impossible, but they 
found that both injection and suction are able to break up streamlines. The total fluid 


mass is also not conserved for the flows of films and drops over porous substrates (Davis & 
Hocking 2000); in fact both drops and films are drawn entirely into the substrate in finite 
time as would be expected. In the analysis of drop evaporation, a number of studies, e.g. 
Anderson & Davis (19951 and more recently Todorova et al. (2012) among others, have 


considered the steady state obtained by imposing injection through the substrate that 
exactly balances the mass lost to evaporation. In the latter study the injection profile was 
imposed according to a Gaussian distribution, and the drop shape is largely independent 
of this distribution as long as the injection is not too large near the drop contact line (note 


that a precursor film model was used by Todorova et al. 2012). For continuous falling 


liquid films over porous substrates there have been linear stability studies invoking a 
Darcy law in the porous medium and a Beavers-Joseph boundary condition at the liquid 
substrate boundary ( Sadiq fe Usha|[2008 Usha et al. 2011). It is found that an effective 
slippage takes place that enhances the instability in the sense that it reduces the critical 


Reynolds number. Slippage models were investigated further by Samanta et al. (2011) 


and an alternative porous medium model is proposed and analysed by Samanta et a,. 


(2013). 


In this paper, we impose blowing or suction through the wall, and perpendicular to it, 
of fluid that is identical to that of the liquid film. The magnitude of the blowing/suction 
is assumed to vary spatially along the planar substrate and hence modifies the no pene¬ 
tration boundary condition there, but we assume that there is no slip along the substrate. 
We envision that such a model would be appropriate to experimental setups where tiny 
slits on the substrate would provide the conduit for fluid to enter and leave the wall. Such 
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Figure 1: Sketch of flow domain showing coordinate system. We consider a fluid layer, 
with mean height h s , bounded along y = 0 by a rigid wall inclined at angle 6 to the 
horizontal, and at y = h(x , t) by a free surface. Fluid is injected through the wall, and 
so the normal velocity at the wall is given by the prescribed function v = F(x,t). 


mechanisms affect the total mass in the film and on physical grounds we can anticipate 
that a net suction would dry the substrate in finite time whereas a net blowing increase 
the total mass and hence the mean thickness at any given time. An increase in thickness 
would consequently increase the local (in time) Reynolds number since it is proportional 
to the mean thickness, hence the flow is expected to become more unstable. In this study 
we will consider the case of blowing/suction that conserves the total film mass (e.g. 
spatially periodic blowing/suction of zero mean), which is possibly the most interesting 
case since it sits on the boundary of the net mass decrease or increase, and hence both 
stabilising and destabilising phenomena can occur depending on the parameters, as will 
be seen later. 

The rest of the paper is organised as follows. In )j2j we discuss the governing equations 
and dimensionless parameters, the scaling and statement of the two long wave models, 
and the choice of the blowing and suction function. The numerical methods used to solve 
these models are described in (|3] In S|4j we explore the steady states and bifurcations 
obtained for non-zero imposed suction, discovering a non-trivial bifurcation structure 
even at zero Reynolds number. We also discuss the distinctive effect of the suction on 
flow streamlines. Linear stability of steady solutions is discussed in fjbj with a focus on 
stability to perturbations of arbitrary wavelength, and thus the effect of suction on the 
critical Reynolds number. In (JbJ we investigate the effect of imposing suction on the 
travelling waves which occur above the critical Reynolds number. In (JTJ we review the 
various initial value calculations performed in this paper, and provide further results. 
Finally, we present our conclusions in 


2. Problem formulation 

2.1. Nondimensionalisation and scaling 

We wish to determine the evolution of a falling liquid film, with mean thickness h s , on a 
slope inclined at angle 9 to the horizontal. We model the liquid as a Newtonian fluid of 
constant dynamic viscosity p and density p, and the air as a hydrodynamically passive 
region of constant pressure p a . The coefficient of surface tension across the air/liquid 
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interface is 7 . We assume that the film flow is two-dimensional, with no variations in the 
cross-stream direction. We use coordinates as illustrated in figure [l| x is the down-slope 
coordinate, and y is the coordinate in the direction perpendicular to the slope, so that 
the wall is located at y = 0 and the interface is defined as y = h(x,t). We denote the 
velocity components in the x and y directions as u and v respectively. 

The dimensionless film flow is governed by the Navier-Stokes equations 

R + V,U X “f Vlly) - Px T 2 T U XX T Uyy, (2.l) 


R (v t + uv x + Wy) = Py ~ 2 cot 9 + V XX + Vyy, 
which are coupled to the incompressibility condition 


u x +v y = 0 . 


( 2 . 2 ) 


(2.3) 


Here we have non-dimensionalised the equations using the mean film thickness h s as the 
length scale, the Nusselt surface speed of a flat film U s = pgh 2 sin d/2p, (Nusselt 19161 
as the velocity scale, h s /U s as the time scale, and pU s /h s as the pressure scale. We 


note that this scaling, based on surface speed, is the same scaling as used by Tseluiko 


et al. (2013) to study the influence of wall topography on the stability of flow down an 


inclined plane. We define the Reynolds number R and the capillary number C based on 
the surface speed U s , so that 


R = 


ph s U s 


C = 


pU s 


(2.4) 


P 7 

We will suppose that the imposition of suction boundary conditions at the wall does 
not alter the no-slip condition, but does affect the no-penetration condition, so that the 
complete boundary conditions at the wall, y = 0 , are 


u = 0, v = F(x , t). 


(2.5) 


At the interface, y = h(x , i), the tangential and normal components of the dynamic stress 
balance condition yield 


(W T ^y) (f ^x) d - (Vy V, x ) 0; 


1 


P~ Pa 


(Vy + u x h 2 x - h x (v x + Uy)) = 


( 2 . 6 ) 

(2.7) 


1 + h% x " v 1 “' x '" x 1 “ y " C (1 + hi) 3 / 2 ’ 

respectively. The kinematic boundary condition on the interface can be written in integral 
form as 

h t - F(x,t) +q x = 0, 

where q(x, t) is the stream-wise flow rate, defined as 

ph 

(2.9) 


( 2 . 8 ) 


q{x, t)= / u(x,y,t)dy. 
Jo 


2.2. Long-wave evolution equations 

We now seek solutions with wavelength L> 1 , and we introduce the long-wave parameter 
6 = 1 /L. We derive two first-order long-wave models, based on an asymptotic expansion 
in the long-wave parameter (a Benney-type model, see Benney|196 6) and on a weighted- 
residual method (following the approach of Ruyer-Quil & M anneville| 2000); derivations 
of both models are presented in appendix |A| We assume that cot 6 = 0(1). To retain both 
inertia and surface-tension effects, we additionally assume that R = 0(1) and C = 0{8 2 ). 
We choose the canonical scaling F = 0(8) so that the imposed suction can enter and 
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compete with the perturbed flow and hence facilitate possible instability enhancement 
or reduction. 

The essential task of the derivation is to estimate the flow rate q(x , t) for a non-uniform 
film. In both models, the mass conservation equation is unchanged from (2.81, and so we 
have 

ht — F + q x = 0. (2-10) 

However, the two models differ in their treatment of nonlinearities in the momentum 
equation, and thus yield different equations for q. 

In the Benney equation (see § A.l I, q is slaved to the interface height h , and is given 

by 


q{x,t) = 


2 h 3 


h? 

y 


2 h cot 0 — 


,V XX 

c 


R 


f 8h 6 h x 2 h 4 F\ 


V 15 


- ■ 


( 2 . 11 ) 


The first-order weighted-residual approach (see § A.2) instead yields an evolution equa¬ 
tion for q: 


:Rh 2 q t + q = 


2 h 3 


/i 3 

y 


2 h cot 9 — 


C 


R 


( 18q 2 h x 
l 35 


34 hqq x 
35 


+ 


hqF 


( 2 . 12 ) 


The Benney and weighted-residual equations are identical when R = 0. For R ^ 0, the 
Benney equation has a single degree of freedom h(x, t ), while the weighted-residual model 
has two degrees of freedom, h(x , t) and q(x, t). Asa result, the weighted-residual equations 
can in principle exhibit richer behaviour. However, despite the additional complexity, the 
weighted residual equations in fact support bounded solutions across a greater range of 
parameter space (Scheid et al. 2004). 


2.3. Choice of blowing and suction function 

In time-dependent evolution, the mean layer height is conserved only if the imposed flux 
function F has zero mean, and hence steady states can only exist if F has zero mean. 
Conserved mean layer height is the natural state for numerical calculations in a periodic 
domain, and is sometimes called ‘closed’ conditions, as there is no net flux out of the 
domain. However, in experiments, closed conditions are not easy to implement, and so 
experimental realisations more typically impose the fluid flux at the domain inlet. In this 
case, there is no direct control over the mean layer height. 

In the absence of blowing and suction, both the open and closed systems support 
uniform flow via the Nusselt solution, and thus we obtain the same scaling whether 
based on the mean layer height or mean flux. In order to investigate the influence of 
suction, we can consider steady states where the mean layer height remains fixed for 
closed conditions, or where there is no change to the mean flux for open conditions. 
However, the bifurcation diagrams for fixed mean layer height correspond most naturally 
to statements about time evolution, as the mean layer height does not change with time. 
Most of the results we present are for fixed mean layer height, but we will also present 
some results for fixed flow rate, and we note that there is a mapping between results for 
these two conditions. 

In the rest of this manuscript, we will consider the simplest functional form for F with 
zero mean, i.e. a single harmonic mode: 


F(x) = Acos(mx) = A cos 



(2.13) 


This function F(x) has the symmetry that the transformation A —> —A is recovered by 
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translation in r by a distance L/ 2, and so translationally-invariant solution measures, 
such as the critical Reynolds number for onset of instability, must be even in A. 

The long wave equations (2.10), (2.11) and (2.12) also apply if F is unsteady, so 
long as F varies on dimensionless timescales no shorter than 0(l/<5), or the dimensional 
timescale h s /(5U s ). Time derivatives of the vertical velocity, and hence F, would feature 
in the equations at the next order in S. We are therefore free, at this order, to impose 
a time-dependence on F, or even choose F in response to the film evolution. The latter 
formulation would be particularly useful in feedback control studies. 


3. Numerical methods 

The numerical calculations that we perform are of three types: computation of steady 
periodic solutions and their bifurcation structure, linear stability calculations of such 
steady states to perturbations of arbitrary wavelength, and nonlinear time-evolution via 
initial value problems. We conduct these calculations using the continuation software 
package Auto-07p (Doedel & Qldman||2009 ) and Matlab. 

The first task, of computing steady solutions, and exploring their bifurcation structure, 
was conducted in Auto-07p by formulating the problem as a boundary value problem 
with periodic boundary conditions. The Auto-07p code is spatially adaptive, and un¬ 
likely to return spurious solutions. We are particularly interested in limit point, pitchfork 
and Hopf bifurcations. The first two of these correspond to bifurcations of steady states, 
and so can be detected and tracked using the same formulation as for standard steady 
states. With regard to Hopf bifurcations, here we are concerned with Hopf bifurcations 
with respect to time, whereby a steady state becomes oscillatory in time when subject to 
a perturbation of fixed spatial wavelength. If the perturbation satisfies periodic spatial 
boundary conditions, the instability will in fact be oscillatory in both space and time. 
We used Auto-07p to track individual Hopf bifurcations by manually augmenting the 
steady system with a boundary value problem for the spatially periodic but non-constant 
eigenfunction, with the temporal eigenvalue determined as part of the solution. 

The linear stability calculations were performed in Matlab, and we used a pseudo- 
spectral method for the spatial discretisation. After spatial discretisation, the governing 
partial differential equation becomes a large system of coupled first-order ordinary dif¬ 
ferential equations, which we can write in the general form 


F(u,u) = 0 


(3.1) 


so that steady solutions Uo satisfy F(uo, 0) = 0. In order to determine the linear stability 
of a steady solution, we suppose that 


u(f) = u 0 + evexp(At), eCl. 


We now expand (3.1) for small e, to obtain 


Jv + AMv = 0, 


<9F 

~ di i 


uo ,0 


A/r 9F 

~ du 


(3.2) 


(3.3) 


u 0 ,0 


which is a generalised eigenvalue problem for A and v, where J is the Jacobian matrix and 
M is the mass matrix. As very few points were needed for the spatial discretisation (we 


typically used 99 equally spaced points), we solved the eigenvalue problem (3.3) directly 


in Matlab. We used Floquet multipliers to determine linear stability to perturbations of 
arbitrary wavelength, and so modified the Jacobian matrix to account for these when 
necessary. 

We note that in the Benney equations, the only time derivatives are those of interface 
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height h, while the weighted-residual equations also feature time derivatives of flux q. 
This means that for the same spatial discretisation, there are twice as many eigenmodes 
for the weighted-residual equations as for the Benney calculations. We found that the 
weighted residual calculations were prone to spurious eigenmodes, which we removed by 
careful comparison of the eigenvalue spectrum at different spatial resolutions. We also 
neglected the neutrally-stable eigenmode corresponding to increasing the total volume of 
fluid in the domain, which arises in both sets of equations. 

Time evolution calculations were always performed in a fixed spatially-periodic domain. 
The spatial problem was discretised via a pseudo-spectral method, while time-derivatives 
were handled via a second-order backward finite difference scheme (BDF2). The resulting 
code is fully implicit, and solved via direct Newton iteration. 

The code was verified by comparing the steady solutions and bifurcation structure 
obtained in Matlab to those obtained in Auto-07p. Further validation was obtained by 
comparison of numerical results to analytical solutions for the shape of small-amplitude 
steady states and to analytical results for the linear stability of uniform and small- 
amplitude states. 


4. Bifurcation structure for steady states 

In the absence of blowing or suction, the only spatially-periodic steady state is a 
uniform film. Introducing periodic suction naturally imposes a spatial structure on the 
solutions, and means that any steady states must be non-uniform. When R > 0 the 
solutions and bifurcations differ between the two long-wave models. The weighted residual 
model avoids the blow-up behaviour sometimes exhibited by the Benney equations, and 
more accurately represents the behaviour of the Navier-Stokes equations at moderate 
Reynolds number. We will generally present results for the weighted-residual model when 
the two models differ, but we note that a non-trivial bifurcation structure emerges even 
at zero Reynolds number. 


4.1. Steady solutions at small A 

We begin by considering the effect of small-amplitude forcing, in the form of blowing and 
suction, on the uniform steady state h = 1. We choose F = Acosmx where m = 2n/L, 
and seek a steady solution for h and q when \A\ <C 1 of the form 


h = 1 + A$t{H) + 0{A 2 ), q=- + A${Q) + 0(A 2 ). 


(4.1) 


The mass conservation equation ( 2.10| ) immediately supplies Q = exp (imx)/(iin). How¬ 
ever, the complex constant H differs between the two long-wave models. The Benney 
result, obtained by linearising (2.111, is 


H i 


Benney 


1 + -imR 

O 


2 8 m 

2 im + -m 2 cot 9 - Rm 2 + —- 

3 15 3G 


exp (imx), 


while the weighted residual equation (2.12) gives 

18 

1 + —imR 
35 


Hwr — 


2 im + Tjm 2 cot 9 — —Rm 2 + 

o oO oO 


exp(imir). 


(4.2) 


(4.3) 


The resulting solutions for H are illustrated in figure [ 2 ] for L = 10 and L = 40. 
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L = 10 L = 40 



Figure 2: The effect of small A on steady solutions with F = Acos(2ttx/L) for L = 10 
and L = 40. The O(A) correction IR(H) is shown for the Benney equations (solid line) 
and the weighted-residual equations (dashed line). The dotted line indicates F(x)/A. In 
both cases, C = 0.05, 8 = 7 t/4 and R = 2. 


The two expressions (4.2) and (4.3) are equal only if R = 0; otherwise both the 


magnitude and phase of H differ between the two equations. At L = 10, the predictions 
obtained via the two models are in reasonable agreement, but they are indistinguishable 
when L = 40. Returning to the variables of the long-wave derivation, we set m = 5M 
and C = 6 2 C , and expand for small 6 ; we find that both models yield 


H = 


1 


2 SiM 


cot 8 

nr 


M 2 R \ 

-~ + -jT +0(5), 

12(7 5 


(4.4) 


which is the expected order of agreement given that terms beyond the second order in 6 


were neglected in the derivation of each model. We note from (4.41 that the magnitude 


of H is inversely proportional to M at leading order, so that for fixed A , the maximum 
perturbation to the interface grows linearly with the wavelength L of the imposed blowing 


and suction. The long wave expression (4.4) also reveals a phase shift between h and F, 
which tends to 7 t/ 2 as S —> 0. 

In order to calculate the small A correction to the mean flux analytically, we must 
expand the steady solution h = H(x), q = Q(x) to 0(A 2 ). For steady states, with 
F(x)Acosmx, we integrate the mass conservation equation ( |2.10| ) to yield 


_, . 2 A sin mx 

Q(x) = x +- 

3 TO 


+ Q 2 A 2 + 0(A 4 ) 


(4.5) 


where we have used the fact that the spatially-averaged flux is even in A. We then expand 
the interface height in A as 

H(x) = 1 + A(pi cos mx + p 2 suitox) + A 2 (p 3 cos 2 mx + p^ sin 2 mx) + 0(A 3 ). (4.6) 


We solve the flux equation, either (2.11) or ( |2.12 1 , at O(A) and 0(A 2 ) to determine the 
constants p 3 , p 2 , p 3 , P 4 and Q 2 . The resulting coefficients are somewhat lengthy and so 
we do not list them here. In the long wave limit to = 6 M. C = S 2 C, with 6 <C 1, both 
the Benney and weighted residual equations yield 


1 


— [ Q(x) dx = - + — 

L L ) 3 + 6 2 V 4M 2 


+ 0(«5 2 ) +0(A 4 ), 


(4.7) 


so that small amplitude, long wave blowing and suction always increases the mean down- 
slope flux. 
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Figure 3: Bifurcation structure for steady solutions to the weighted residual equations as 
A increases, subject to fixed mean layer height of 1, and calculated in a periodic domain 
of length 2 L. Here F = Acos(2nx/L), L = 10, C = 0.05, 9 = 7t/ 4, with R = 0,3,6,9 in 
(a-d) respectively. The shaded inset solutions lie on the dashed branch of subharmonic 
steady solutions, which is unstable; all other solutions are harmonic, with period L. 
Solutions filled in black are stable, while white solutions are unstable. We also indicate 
pitchfork bifurcations at A = Ap (♦), limit points at A = App{%), Hopf bifurcations 
within the domain of length L (■), and states with ming = 0 (★). 


4.2. Steady solutions as A increases 

As A increases, the film profile deviates from the uniform state, and the nonlinearities in 
equations (2.11) and (2.12) can lead to bifurcations between steady solutions. Figure [3] 
shows the behaviour of steady solutions to the weighted-residual equations as A is varied 
for a selection of R at fixed L , 9 and C subject to constrained mean layer height 1. 

For each R, the only steady solution at A = 0 is the uniform state with h = 1. All steady 
solutions for A > 0 are non-uniform, and the extent of this non-uniformity increases with 
A when A is moderate. However, for sufficiently large A , there are no steady solutions 
that describe continuous liquid films, due to a limit point in the bifurcation diagram. 
We can follow the steady solution branches around the limit point, and find that each 
branch eventually terminates when the minimum film height tends to zero, so that the 
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Falling liquid films with blowing and suction 


(a) A = 0.8, n = 1 


(6) A = 0.72, n = 2 




Figure 4: Time evolution from a near-flat initial state, for R = 0, L = 10, C = 0.05, 
F = Acos(27ra;/-L), in (a) a domain of length L for A > Ap so that there are no steady 
states, and (6) a domain of length 2 L with A > Ap so that all steady states are unstable. 
In both cases, the film height vanishes in finite time at a single point in the domain. 

layer dries up at some position. The film height profiles for these drying states are shown 
as insets in figure [3] 

When A is sufficiently large that there are no steady solutions, we explore the system 
behaviour by conducting time-dependent simulations in a periodic domain of length L, 
with a typical result shown in figure [4^ a). We find that the film thins rapidly, and is 
able to dry in finite time due to the fixed speed of fluid removal. We have also conducted 
unsteady simulations starting from slightly-perturbed states on the solution branch with 
the lower minimum height. We find that these steady states are unstable, and the in¬ 
stability manifests as thinning and drying behaviour, rather than tending towards the 
linearly stable steady states with a larger minimum film thickness. 

The solutions shown in figure[3]are all computed subject to a fixed mean layer height of 
1, and are plotted according to the minimum value of h in a single period. An alternative 
solution measure is the mean flux q, averaged over one spatial period. Figure [5^ a, c) 
shows the bifurcation structure plotted in terms of the mean flux; we find that imposing 
blowing and suction at finite wavelength can either increase or decrease the mean flux. 
Figure[5](&, d) shows bifurcation diagrams for two values of R with fixed q = 2/3 and mean 
height free to vary (this condition is appropriate to experimental investigations performed 
under open conditions, corresponding to fixed mean down-slope flux). We find that, as 
was the case for the fixed h = 1 calculations, there is a limit point corresponding to a 
maximum value of A with steady solutions. However, time-dependent calculations for 
open conditions require non-periodic boundary conditions, which we have not pursued 
here. 


4.3. Subharmonic steady states 


The steady solutions discussed in § |4.2| are limited to cases where both the steady so¬ 
lution and the imposed blowing and suction are periodic with the same wavelength L. 
However, subharmonic steady solutions may also exist, for which the steady solution is 
spatially periodic with period nL : where n is an integer equal to 2 or greater; we calculate 
subharmonic steady states by continuation in a domain of length nL. For each case in 
figure [3j and also for the fixed q = 2/3 calculations shown in figure [5] we have detected 
subharmonic solution branches with n = 2. In each of these cases, the subharmonic so¬ 
lutions emerge via a subcritical pitchfork bifurcation at A = Ap in the steady solution 















Figure 5: The bifurcation structure for steady solutions to the weighted residual equations 
as A increases, subject to fixed mean height (a, c) and fixed mean flux ( b, d). Here L = 10, 
C = 0.05 and 9 = 7r/4. The symbols mark the same bifurcations as in figure [3j but we 
do not show Hopf bifurcations here. Solution branches terminate at a point where the 
minimum layer height vanishes. 


structure, and the subharmonic steady states are all unstable. The unstable eigenmode 
of these states results in a drying process similar to that occurring for unstable harmonic 
steady states. 

As the subharmonic steady states shown in figure [3] are unstable, the only possible 
stable steady states are harmonic, with period L. However, for A > Ap, the period-L 
steady states are also unstable to perturbations of wavelength 2L; this is a consequence of 
the subcritical nature of the pitchfork bifurcation. The subharmonic instability eventually 
results in a film thinning and drying event, but at only one of the local minima within 
the domain of length 2 L. One such drying sequence is illustrated in figure |4^&). 


4.4. Streamfunction and flow reversal 

In addition to altering the interface height, the imposed suction also fundamentally affects 
the arrangement of streamlines within the fluid film. As the flow is incompressible, we 
can write the velocity in terms of a streamfunction ip(x , y) such that (it, v ) = (ip y , —ipx)- 
Under the weighted-residual formulation, we have 

u= ^ = f{h-h) +oiS) - v= (48) 


with boundary condition v = F(x) on y = 0. Recalling that F(x) = q'(x) for steady 
solutions, we can write the solution for ip, up to the addition of an arbitrary constant, as 


iP 


-q(x) 


3 q(x) f 
h{x) \2h 



+ 0 ( 6 ). 


The same result applies to 0(6) in the Benney equations. 


(4.9) 
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According to (4.9), there are no stagnation points in 0 < y < h if q > 0, and points 


along the wall are stagnation points if and only if F(x) = 0. In the absence of suction, 
every point on the wall is a stagnation point. With non-zero blowing and suction, there 
are isolated stagnation points on the wall at the zeros of F(x). 

We now examine the steady flow near such an isolated stagnation point, located at 
x = Xp, y = 0, with F(x o) = 0. As q'{x) = F(x ), we can write 


q ~ q(x 0 ) + F\x o) — ^ + 0{x - z 0 ) 3 . 


(4.10) 


Substituting (4.10) into (4.9) yields 


if - q{ x 0 ) - 


F'{x 0 )x 2 


3q{xo)y 2 n(( , 2 


-( y-vo ) 2 )- (4.11) 


If F’ (xp)q(xp) > 0, the stationary point at Xp is a saddle point of if, and the streamlines 
are locally straight lines, such that 


y 2 F' {xp)h{xpf 


3 q{xp) 


(4.12) 


The streamlines become increasingly vertical as q(Xp) —> 0. In contrast, if F' (xp)q(xp) < 
0, the stationary point at xp is a local extremum of if, and the streamlines are locally 
elliptical, with 


- F'(xp)x 2 + 


3q(xo)y 2 

h(x 0 ) 2 


= const. 


(4.13) 


For any smooth, periodic, non-zero F(x) with mean zero, there must be some stagnation 
points Xp with F'(x o) > 0 and others with F'(xp) < 0. There are no steady solutions in 
figure [3] with negative h, but the same is not true for q; for each R with solutions shown 
in figure [3j there is an amplitude A* above which all steady solution have q(x) < 0 over 
a finite interval in x. As q x = F, the minimum of q occurs at a point xp with F(xp) = 0 
and F'(xp) > 0. 

There are no stagnation points inside the fluid domain, and so streamlines never cross. 
For small A, q(x) > 0 for all x, and so stagnation points with F'(x) > 0 are saddle points 
of if, while those with F'(x) < 0 are extrema of if. Two examples of the flow field can be 
seen in figure [b](a, b). Streamlines emanate into the fluid domain from each stagnation 
point with F'(x) > 0, and these stagnation points are connected by a streamline which 
separates the fluid into a layer in which fluid particles propagate down the plane, and a 
layer in which fluid particles must enter and leave the flow domain via injection through 
the walls. At A = 0, all particles propagate. As A increases, the propagating layer thins 
and the injection layer thickens. 

At A = A*, the stagnation point at x = xp, y = 0 switches from a saddle point to an 
extremum of if. The corresponding change in the flow held is illustrated in the transition 
from figure [6](&) to[6]) c). For A > A*, q(x) is negative for some x, and hence by ( |4.8[ ) the 
horizontal velocity is directed up-slope for all y. As a result, no particles can propagate 
more than one period down the plane, and so the propagating layer vanishes for A > A*. 
The width of the region with up-slope how increases with A, until eventually there are no 
further steady solutions. Nonlinear time evolution calculations at large A show that the 
him dries at a point just upstream of the negative-g stagnation point. An instantaneous 
snapshot of the system at the moment of drying is shown in figure [6](e). 
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(b) Steady state for A = 0.4 



(c) Steady state for A = 0.5 



(e) There are no steady states for A = 0.8; this is a final snapshot before drying. 


Figure 6: Solutions for the interface shape, velocity held and stagnation points for 6 = 
7t/4, C = 0.05, L = 10, R = 0. Periodicity is enforced with L = 10, but solutions are 
plotted over two periods, and are shown with aspect ratio 1. Stagnation points with 
q(x)F'(x) < 0 and with q(x)F'(x) > 0 are indicated by O and □ respectively. In 
( a-b ), q > 0 everywhere, so a streamline emanating from the stagnation point with 
q(x)F'(x) > 0 divides fluid into particles which never meet the wall (yellow), and particles 
which enters and leaves through the wall (blue). For A > A* = 0.46, all steady solutions 
have regions of negative q , corresponding to a region of upstream how near the stagnation 
point (between the vertical lines), and all fluid particles must reach the wall. 


5. Linear stability 

5.1. Stability of uniform state 

In the absence of imposed suction, the only steady state with mean layer height unity 
is the Nusselt him solution, with h = 1 and q = 2/3. We linearise about this state, and 
seek eigenmodes proportional to exp (ikx + At). The Benney model yields the eigenvalue 
A directly, as 

k 4 
‘ 3 C' 


, , 8 R 2 

A = —2 ik + At (-cot t 

15 3 


(5.1) 


For the weighted residual model, the linear stability analysis yields a quadratic equation 
for A: 


— A 2 + A 
5 


„ , 68 R\ l2 

1 - ik—— + 2 ik + k 2 
105 / 


2cotd 


^-^ 1 = 0 . 

3 C 35 


(5.2) 
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In both models, the stability threshold for fixed k is 

R < R h = ~ cot 9 + J^zk 2 . 
4 8 C 


(5.3) 


The linear stability threshold (5.3) is also valid for perturbations with long wavelengths 
in full solutions to the Navier-Stokes equations (Benjamin 1957). At R = Rh, both (5.1) 
and (5.2) have a root with negative real part for R < Rh, the value A = Ao = —2 ik at 
R = Rh, and positive real part for R > Rh- We note that (5.2) also has a second root 
for A, but this root always has negative real part when R > 0. 


5.2. Perturbations of arbitrary wavelength about a periodic base state 

When the base state is spatially uniform, the eigenmodes can be written as 3?(exp(ifccc + 
At)), and by calculating A for all real k, we can determine linear stability to perturba¬ 
tions of all wavelengths. However, when non-zero suction and blowing is imposed, the 
base state is not uniform, and so the eigenmodes are no longer simple exponential func¬ 
tions. Fortunately, Floquet-Bloch theory allows us to compactly describe eigenmodes of 
arbitrary wavelength. 

We suppose that we are given a steady solution h, = H(x), q = Q(x) to the forced 
equations, and consider the evolution of arbitrary small perturbations, so that 

h(x,t) = H(x) + eiR{h(x,t)} + 0(e 2 ), q(x,t) = Q{x) + elR{q(x,t)} + 0(e 2 ), (5.4) 


with e<l. The mass conservation equation ( |2.10 1 becomes 

ht + q x = 0. 


At 0(e), the Benney flux equation (2.11) yields 

q = b 0 {x)h + bi(x)h x + b 2 (x)h xxx 


(5.5) 

(5.6) 


where 


b 0 = H 2 [ 2 — 2H X cot 6 + 


H x 


C 


^H x - —H 4 F, 
5 3 


(5.7) 


h = -~^F[ 3 cot8+ — RFl 6 , b -2 = 
3 15 


H 3 
3O' 


(5.8) 


The weighted residual equivalent of (5.6) additionally involves a term proportional to q t - 
Given a base solution H(x), Q{x) and F{x), the coefficients 6o, &i an d b 2 are all known 
periodic functions of x, with the same period as the base solution. 

We now invoke the Floquet-Bloch form, observing that as (5.5) and (5.6) are linear 
equations with coefficients that are periodic in x with period L 1 the eigenfunctions can 
be written as 


h(x,t) = e M+lkx h k {x), q{x,t) = e xt+lkx q k (x), (5.9) 

where the Floquet wavenumber k is real, and hj- and qk are periodic functions of x with 
period L. Setting k = 0 recovers perturbations of wavelength L, while very small but 
non-zero k corresponds to very long wave perturbations. The solution is stable to pertur¬ 
bations with period L if 3?(A) < 0 for all eigenfunctions when k = 0. The base solution 
is linearly stable to perturbations of all wavelengths if 3?(A) < 0 for all eigenfunctions for 
each real k £ [0,7r/L], 


5.3. Effect of small-amplitude blowing and suction on stability 
We now analyse the effect of small amplitude blowing and suction in the form F = 
Acosmx on the stability of eigenmodes with underlying Floquet wavenumber k, which 
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will allow us to determine how such forcing affects the critical Reynolds number. We do 
so by expanding the eigenvalue A for R close to Rh, and for small A: 

d\ 

A = A 0 + {R-R h ) ZA 2 + 0((R - R H f) + 0((R - R H )A 2 ) + 0(A 4 ). (5.10) 

Oix 

We are concerned with eigenvalue equal to Ao = —2 ik at R = Rh, which is a root of 
both the Benney and weighted residual characteristic equations. We can evaluate the 


term dX/dR which appears in (5.10) by differentiating (5.1) or (5.2) as appropriate. 
The eigenvalue is even in A because the transformation A —> —A can be recovered by 
translation in i by a distance L/2, and the original eigenmode exp (ikx + Xt) has no 
preferred position; thus the leading order contribution for small A is 0(A 2 ). 

The effect of the imposed suction is encapsulated in the coefficient Z, which depends 
on the suction wavenumber to, the perturbation wavenumber k , and C and 6 . In order to 
calculate W and Z, we need to calculate both the base solution and the eigenfunction to 
0{A 2 ). The calculation of the base solution to 0(A 2 ) was discussed earlier in the context 


of steady states, and the required expansion is given by (4.51 and (4.6). 


For small-amplitude steady solutions, the suction function F and the base solution H , 
Q are all periodic with wavenumber to, and so the coefficients in the linearised equations 


(5.5) and (5.6) are also periodic with wavenumber m. As a result, all eigenfunctions of 


(5.5) and (5.6 1 can be written in the Floquet form given by (5.9). The unknown functions 


hk and q k are periodic with wavenumber to and are constant when A = 0, and can be 
expanded for small A as 


h k (x) = 1 + A(C X e imx + C 2 + C 3 e~ imx ) 

+ A 2 (D 1 e 2imx + D 2 e imx +D 3 + D 4 e~ imx + D 5 e 


— 2 imx\ 


and 


q k {x) = ^ + A(E ie imx +F 2 + E 3 e~ imx ) 
k 

+ A 2 (G 3 e 2imx + G 2 e imx + G 3 + G 4 e~ imx + G 5 e~ 2imx ), 


(5.11) 


(5.12) 


where the complex constants G\, C 2 , C 3 , D\, D 2 , D 3 , ZAt, D 5 , E\, E 2l E 3 , G i, G 2 , G 3 , 
G\ and G 3 are to be found along with Z. We must also choose a normalisation condition 
on the amplitude and phase of the eigenvector, and for this we use the condition 


MO) = i. 


(5.13) 


To determine the unknown constants, we set R = Rh, solve for the steady state given by 

(4.5) and (4.6), and then substitute the eigenvalue expansion from ( 5.10| ) and the eigen¬ 
function from ( |5.9[ ), ( |5.11 ), (5.12) into the two linearised equations (5.51 and (5.6), and 
also the normalisation condition (5.13). We solve the linearised equations first at O(A) 
and then at 0{A 2 ), at each order obtaining linear systems for the unknown coefficients 
including Z. We perform this calculation in Maple, and obtain a lengthy expression for 
Z as a function of to, n , C and 9 ; the full result depends on whether we have used 
the weighted-residual or Benney equations. Once Z is known, we can obtain the neutral 
stability curve for fixed k , defined by the condition 5i{A} = 0, so that 


R~R h + A 2 W, W = -3f£{Z} 


\dR 


-1 


(5.14) 


Figure [T] shows the effect of small A on the stability boundary, as quantified by the 
weighted-residual results for W. For 9 = 7r/4, we find that forcing via blowing and suction 
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Figure 7: Effect of small A blowing and suction with wavelength L on the critical Reynolds 
number for perturbations of underlying wavenumber k as measured through the weighted- 
residual results for W, with C = 0.05. The filled circles on each curve mark the point 
where k = 2i t/L. The dashed lines show the long-wave asymptotic estimate (5.16|. 


at wavelength L = 2tt/ m = 10 destabilises short wavelength perturbations, but has little 
effect on long wavelength perturbations. However if the forcing wavelength L > 20, W is 
positive for all k, and so forcing has a stabilising effect on all perturbations. Furthermore, 
for long wavelength forcing, the magnitude of W increases with L , and the minimum value 
of W occurs at k = 0; this value is positive and increases rapidly with L. For 0 = 7t/ 2, 
forcing with wavelength L has a destabilising effect on short waves for L = 10, and a 
stabilising effect for L = 20, 40, 80. However, W tends to approximately —6 as k —> 0, and 
so the imposed forcing in fact slightly destabilises long wave perturbations. For 9 = 37t/ 4, 
the behaviour for small perturbation wavenumber is reversed from the 9 = 7t/ 4 case; long 
wave forcing is always destabilising for k = 0, and becomes increasingly destabilising as 
L is increased. 

In figure |8j we survey the sign of W for k = 0 and k = m, with m and C varying, for 
three values of 9 , for both the Benney and weighted-residual models. The two models 
concur that for each inclination angle 9 , there are some regions of (C, to) parameter 
space where the imposition of blowing and suction stabilises the uniform state, and other 
regions where such forcing destabilises the flow. For 9 = 7t/ 4, both models predict that 
long wave forcing, with to —> 0, is stabilising to perturbations with Floquet wavenumber 
k = 0 and with wavenumber k = in. For 9 = 7t/ 2, forcing is destabilising in the limit 
k —> 0 for all C and to, but has a stabilising effect on perturbations with k = m if m is 
not too large. For 9 = 37r/4, the long wave forcing is destabilising to perturbations both 
with k = 0 and with k = m. 

The critical Reynolds number increases quadratically with perturbation wavenumber 
k in the absence of blowing and suction, and so for very small A, the lowest critical 
Reynolds number must also be obtained at small k. We now consider W in the limit of 
long-wavelength blowing and suction and long perturbation wavelength, setting m = SM , 
k = 6K and C = S 2 C , and expanding the full expression for Z with 6 <C 1. Both the 
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Benney, k = 0 WR, k = 0 Benney, k = m WR, k = in 



Figure 8: The effect of small amplitude blowing and suction with wavenumber m on 
perturbations with k = 0 and with k = m in the Benney and weighted-residual (WR) 
models, as described by the quantity W defined in ( 5.14| ); W > 0 in the shaded regions, 
so small-amplitude blowing and suction increases the critical Reynolds number, and has 
a stabilising effect on the flow. Results for 6 = 7r/4, 7 t/2, 37t/ 4 are shown in rows (a-c) 
respectively. 


Benney and weighted-residual results yield 




3 iK 
ASM 2 


K 2 cot 9 
2 M 2 


K 2 
6 C 


4 1 


11 K 
12CM 2 


+ 0(S 2 ), K 


fdX 


8 K 2 S 2 

15 


Substituting (5.15) into (5.14) yields the leading order long-wave expansion 

11 A' 2 


0(S 4 ). 

(5.15) 

- 0 ( 1 ). 

(5.16) 

The prediction (5.16) is plotted in figure [7] for comparison to the non-long-wave results 
obtained by direct evaluation of the Maple expressions for Z and dX/dR. 

For small-amplitude blowing and suction at fixed wavenumber m , the critical Reynolds 
number for instability to perturbations of all real wavenumbers k is 


W = - — 
8<5 2 


cot 9 
2M 2 


1 

60 


12 CM 2 


o(i)=- 1 4 


cot 9 | 1 life 2 

2^2 + 6C~ 12 Cm 2 


R*(m) = min 
fee R 


5 „ 5k 2 

r ote+ w 


A 2 W(m, k, C, 9) 


Using the long wave expansion for W given by (5.16), we obtain 

5 A 2 


5 15A 2 

R* (to) = - cot 9 + ——- cot 9 - —— + mm 
4 16 to 2 160 fees 


k 


55^4 2 \ 


80 + 32to 2 0 J 


(5.17) 


(5.18) 


The term in square brackets is positive, and so the minimum of the expression occurs 




































































Falling liquid films with blowing and suction 


19 


at k = 0. As expected, this implies that long-wave perturbations are the first to become 
unstable as R is increased, both for A = 0 and also for small A in the long wave limit. 
Proceeding with the long wave approximation, we evaluate the minimum of ( 5.18| ) as 

< 5 ' 19 ) 


R* 


cot 9 


15 


16 V m 2 


while if perturbations restricted are restricted to those with wavenumber k = m, we find 

/ cot 6 3 


R 


m—k 


5 „ 5m 2 


_ 

16 I m2 


2 C 


(5.20) 


The stability boundaries for perturbations with all k and with k = m are plotted in figure 
[I for nonlinear solutions to the weighted-residual equations, small A predictions obtained 
using the full version of W, and the long-wave, small A predictions (5.19) and (5.20). 


We see that the first two boundaries are in good agreement with each other when A is 
sufficiently small, for both L = 40 and L = 10. The long-wave small-amplitude predictions 
are in good agreement with the other two versions of boundaries when L = 40, but are 
poor agreement when L = 10. 

In the case 9 = n /4, R* is positive in the absence of suction, and long-wave blowing 
and suction can either decrease or increase R*, depending on the values of 9, m and C, 
and here it is reasonable to consider an optimisation strategy. To maximise the stabilising 


effect at fixed A for 9 < 7t/ 2, we should choose m as small as possible, and in fact (5.19) 
predicts that R 


+00 as m 


0 with A fixed. Regarding constraints on the magnitude 
of film height variations, we note that these are of order A/m for A <C 1 (see § 4.1) , and 


so we can approximately constrain height variations by keeping a = A/m fixed. We find 
that the largest R* is again obtained as m —> 0, but in contrast to the fixed A case, R* 
is bounded as to 0. 

The governing equations are valid for 0 < 9 < n, with 9 = 7r/2 corresponding to 
flow down a vertical wall, and w/2 < 9 < n corresponding to flow along the underside 
of an inclined plane. In the absence of blowing and suction, R* = 0 for 9 = 7r/2, and 
R* < 0 for 7 t/2 < 9 < tt. Negative Reynolds numbers are not physically relevant, but 


consideration of (5.19) may still give some indication of the effect of blowing and suction 


when 9 > 7r/2. According to (5.19), when 9 > 7r/2, imposing blowing and suction at long 


wavelengths always reduces f?*, and so is destabilising. Our other results for 9 = 37t/ 4, 
shown in figures [7] and [8j also suggest that suction destabilises long wave perturbations. 


5.4. Finite A stability results 

Figure [9] shows stability regions as R and A vary, at fixed values of C and 9 , for imposed 
suction with wavelength L = 10 and L = 40. For L = 10, introducing finite amplitude 
blowing and suction decreases the critical R for linear stability to perturbations of wave¬ 
length L and to perturbations of arbitrary wavelength. In contrast, when the imposed 
suction has wavelength L = 40, increasing A from zero initially increases both of these 
critical Reynolds numbers, and so has a stabilising effect on the base flow. 

The boundary of the region stable to perturbations of fixed wavelength L is defined 
by a Hopf bifurcation. When A = 0, this is the Hopf bifurcation at R = Rh , for which 
the eigenmode can be written as h = e lkx with k = 2-n/L. We can use Auto-07p to 
track the Hopf bifurcation as A varies. As A is increased from zero, the eigenmode is 
modulated and is no longer a single Fourier mode, but remains periodic with period L. 

Calculation of the stability boundary in the presence of patterned blowing and suction 
for perturbations of all wavelengths is not as simple as the case of tracking a single 
bifurcation, as we do not know a priori which wavelengths are most unstable. However, 
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(a) L = 10, with C = 0.05, 9 = 7r/4 



(b) L = 40, with C = 0.05, 9 = 7t/4 


Figure 9: Stability properties for steady solutions as R and A vary at fixed 9 = 7r/4 
and C = 0.05, for the weighted-residual model. Steady solutions are divided into three 
linear stability categories as indicated by the legend, with finite-H stability boundaries 
shown by solid black lines. Two small-H estimates for these boundaries are shown: the 
full result correct to 0(A 2 ) (dash-dotted line), and the long-wave, small A estimates 
(5.19) and ( 5.20| ) (solid white line). We also indicate solutions with min(g) = 0 (dashed 
line), and the pitchfork bifurcation Ap (dotted line). 


in the absence of blowing and suction, we know that long waves, with 0 < k 2 <C 1, are the 
first to become unstable as R is increased. For R > Rh, there is a wavenumber cutoff, 
so that perturbations with k 2 < k 2 are unstable, and there is a finite k 2 with maximum 
growth rate 3?(A). However, it is possible that for larger A, imposing blowing and suction 
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Figure 10: Finite A stability results across a range of blowing and suction wavenumber 
m, with F = Acosmx, C = 0.05 and 9 = 7t/ 4. Each black indicates that for the given Ft 
and m, there is some A for which there is a steady state stable to linear perturbations of 
all wavelengths. The dashed line is obtained by explicit tracking of the maximum stable 
R for fixed Floquet number k = 0.01; this curve follows the stability boundary relatively 
well until it diverges at to « 0.05. Along this tracked curve, a = A/m is between 0.5 and 
1. The blue line shows the small-A, long-wave stability prediction (5.191 with a = 1, but 
this prediction underestimates the range of stable R for small to. 


can destabilise at some wavelengths while stabilising at others, and so the system may 
first become unstable at a non-zero wavenumber. For finite A we conduct a brute force 
computation of stability properties over a large number of perturbation wavelengths to 
determine stability with respect to perturbations of all wavelengths. 

For the cases shown in figure [9] we find that at small A , the boundary for stability to 
perturbations of arbitrary wavelengths agrees well with the asymptotic results derived in 
the previous subsection. For both L = 10 and L = 40, this boundary connects smoothly 
to a finite suction amplitude A at R = 0. There is a turning point with respect to R in the 
L = 40 results, and so the largest R at which the system is stable occurs for a non-zero 
A. However, for L = 10, there is no turning point, but there is a second ‘island’ region 
where the steady state is stable to perturbations of all wavelengths. This region requires 
moderately large A and R. Results from time dependent simulations in and around the 
island are shown in figure [14] As the island region occurs only for L = 10, it is possible 
that it is simply an artefact of applying long-wave models at relatively short wavelengths. 

5.5. Optimal wavelength for blowing and suction to obtain a stable steady state 

Determination of the largest R that can be stabilised in an infinite domain by imposing 
a suitable suction profile is a complicated question, requiring finite amplitude results 
for the steady solutions and for their linear stability operator. We can address this task 
numerically, by choosing a suction wavenumber to, increasing R slightly from Rq = 
1.25 cot 0, and testing the stability of steady solutions beginning from A = 0 until A is 
too large for steady states. Figure [To] shows numerical results from this search for the 
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case 8 = 7t/4, C = 0.05, using the weighted-residual model. We find two ranges of m in 
which imposing suction can increase the critical Reynolds number. 

Firstly, at large wavelengths with L ss 200, imposing blowing and suction allows the 
critical Reynolds number for stability to all perturbations to increase from the unforced 
value of 1.25 to around 8. However, figure [To| seems to indicate a downturn in the max¬ 


imum stable R at very small m. From the small A, long-wave result (5.19), we expect 


forcing via blowing and suction at small m to have some stabilising effect, but this result 
does not predict the magnitude of the stabilising effect on its own, as we do not know A. 
Furthermore, estimates of the maximum stable R typically require large A calculations, 
and so predictions based on the long-wave, small-amplitude expression (5.19) are not 
particularly useful in this case. 

The second region in which steady solutions are stable to perturbations of all wave¬ 
lengths in figure [lO] appears when the wavenumber m for blowing and suction is relatively 
large, and only at moderately large R. This region corresponds to the island region in 
the L = 10 results in figure [9j but there is no such stable island in the results for L = 40. 


6. Travelling waves 

6.1. Travelling waves in the absence of suction 

In the absence of suction, the system can support travelling waves, propagating at a 
constant speed U without changing form. We can write the solution as 

h(x,t)=H{ C), q{x,t) = Q{ C), C = x-Ut, (6.1) 

where U is the unknown constant wave speed. H and Q must satisfy 

- UH' + Q' = 0 (6.2) 


and 

- UQ' + Q = ^ ^2 - 2 H' cot 8 + + R Q 2 H' - , 


(6.3) 


where a prime indicates a derivative with respect to C,. If H((f) is periodic with period 
L , then the travelling wave solution is spatially periodic with period L , and temporally 
periodic with period T = L/U. We can compute large amplitude travelling waves nu¬ 
merically; figure 11 shows one such travelling wave solution for 8 = 7r/4, C = 0.05, 
R = 3. 


Individual travelling waves may be stable or unstable, and branches of travelling waves 
can undergo a range of bifurcations. However, small amplitude travelling waves connect 
to the uniform film state via a Hopf bifurcation at R = Rh- This Hopf bifurcation can be 
supercritical or subcritical, with stable, small-amplitude travelling waves observed near 
the bifurcation only in the supercritical case. We can determine the criticality of the 
Hopf-bifurcation by solving for small-amplitude limit cycles near to the critical Reynolds 
number via the following expansion : 


C = x — Ut , H(() = 1 + ecoskif + e 2 H 2 (() + 0(e 3 ), (6.4) 


i?2(C) = r i cos 2&C + r 2 sin 2fc£, q{x,t) = UH(() + S, (6-5) 

U = U 0 + e 2 U 2 + O(e A ), S = S 0 + e 2 S 2 + 0(e 4 ), R = R H + e 2 R +0(e 4 ). (6.6) 

We expand the equations for e <C 1, and must solve the equations at up to 0(e 3 ) in order 
to determine R. We find that small amplitude travelling waves always travel downstream, 





23 


Falling liquid films with blowing and suction 


with speed U = Uq = 2, which is twice the velocity of particles on the surface. The 
bifurcation is supercritical if R > 0 and subcritical if R < 0. The Benney equations yield 


- 120C 2 - 60 C 2 k 2 cot 2 9 - 120 Ck 4 cot 9 - 45 k 6 

R ~ 64Cfc 4 


(6.7) 


which may be positive or negative. However, for the weighted-residual equations, we find 


- _ 4410C 2 + 6670fc 2 C 2 cot 2 9 + 12235fc 4 C7cot 9 + 4450/c 6 
~ 2352 Ck 4 


( 6 . 8 ) 


When 9 < 7t/ 2, ( |6.8[ > yields R > 0 and so the Hopf bifurcation in the weighted residual 
model is supercritical. In the long-wave limit, with k = 8K and C = 8 2 C , both (6.71 and 


(6.8) yield 


R 


15 C 
8<5 2 A' 4 


+ 0(S 2 ) = ^ + °(S 2 ) 


(6.9) 


which is positive. 

The quantity R determines the criticality of the Hopf bifurcation, and so also governs 
the stability of small-amplitude travelling waves in the neighbourhood of the bifurcation. 


travelling waves. 


However, the branch of travelling waves may undergo further bifurcations (Scheid et al. 


2004), and so the value of R does not necessarily prescribe the stability of finite-amplitude 


6.2. Influence of heterogeneous blowing and suction on travelling waves 
Introducing spatially-periodic suction means that the system is no longer translationally 
invariant, and so we cannot obtain true travelling waves for non-zero A. For the final part 
of our analysis, we calculate the effect of small amplitude suction on travelling waves, and 
consider in particular how travelling waves can transition to steady, stable, non-uniform 
states as the amplitude of the imposed blowing and suction is increased. 

We now perturb the travelling wave by introducing a small-amplitude suction F = 
Af{x) with f(x) periodic with period Lp, and A small. We expect to find 

h(x,t) = H(() + Ah(x,t) + 0(A 2 ), q{x,t) = Q(() + Aq(x,t) + 0(A 2 ), (6.10) 


where the behaviour of h and q is in some way related to the periodicity of the original 
travelling wave with respect to f with period L, and of the blowing and suction function 
with respect to x with period Lp. 

The weighted-residual equations for h and q yield, at O(A), 

h t -f(x) + q x = 0 (6-11) 


and 

q + w 6 (()q t = w 0 (()h + WxiOhx + w 2 (()h xxx + w 3 (()q + w 4 (()q x + w 5 (()f(x), (6.12) 


where 


w 0 = H 2 (2- 2 H' cot 9 + 


H" 


34 RQQ' 


Wi = 


2 rr3 n 18I?Q 2 
--H 3 cot9+— 

3 35 


C J 35 

/ 36QH 1 34 HQ' 


H 3 

W2 = 3C’ W ’ = R \ 35 


35 


34 RHQ RHQ 2 RH 2 

m = -35-, «* = —, -6 = — 


(6.13) 

(6.14) 

(6.15) 


The coefficients Wi,i = 0, ..., 6 are functions of f. If we set f(x) = 0 in these equations, we 
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recover the equations governing the evolution of linearised perturbations to the travelling 
wave, i.e. the equations of linear stability. For non-zero /, we can integrate forward in 
time from a given initial condition h(x, 0) = ho(x), q(x, 0) = qo{x), to determine h(x,t) 
and q(x,t). 


The system (6.11) and (6.12) features some terms that are known functions of £, and 


others that are known functions of x , but the system is autonomous with respect to t. 

terms of 

qx + q(- 


We therefore choose to rewrite the equations (6.11) and (6.12) in terms of x and £: 

(6.16) 


ht—Uh^, h x -^h x + h^, qt^f—Uq^, q x 

We obtain a system in two variables, x and £, with equations 

-Uh{-f(x) + q x + q{ = 0 (6-17) 

and 

q-Uw 6 (C)qc = w 0 (C)h + wi(C){h x + h (: ) + W2{0(h xxx + 3h xxt ; + 3h xii: + h(; (: (;) 
+w 3 (( )q + w 4 (C )(Qx + Qc) + w 5 (C )f(x). 

(6-18) 

We now regard £ and x as independent variables. Under this transformation, x remains 
as a purely spatial variable, but ( has a dual identity, incorporating both spatial and 
timelike components. The statement of initial conditions becomes more complicated in 
the new variables: 


h(C = x) = h 0 (x), q(( = x)=q 0 (x) 

so that the initial conditions are spread across the whole range of (. 
If the Fourier expansion of f(x) is 

f(x) = f m cos mx + g m sin mx, 


(6.19) 


( 6 . 20 ) 


we can write the general solution of (6.17) and (6.18) as 


h = ^ Jm{0 cos mx + K m ((') sinma; + h*((, x), (6-21) 

m 

q = ^ M m (Q cos mx + N m (() sin mx + q*{x) (6.22) 

m 

where the functions J m ((), M m (Q and N m (£) are periodic in (. The remaining 

terms, h* and q*, satisfy the homogeneous versions of (6.17) and (6.18) obtained by 


setting f(x) = 0, which are exactly the equations for linear stability of the underlying 
travelling wave. We can calculate the periodic functions J m , K m , ..., corresponding to a 
limit cycle, for any periodic travelling wave. However, we would only expect to observe 
the limit cycle behaviour in initial value calculations if the limit cycle is stable. If the 
travelling wave is stable, all solutions h*, q* to the homogeneous problem eventually 
decay to zero and so h and q are limit cycles, periodic in time. 

We now calculate h for a small-amplitude limit cycle driven by f(x) = cos mx , so that 


the forcing Fourier series has only a single component. The sums in (6.21) are then over 
a single value of to, for which we obtain 

— U J' + mN m + M' = 1, 


- UK' - mM m + NL = 0, 


(6.23) 

(6.24) 











M m - w 6 UM' m 


and 

Nm - WQUN' m 
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w 0 J m + Wi{j'm + fnKm) + w 2 (J"' + 3 toA" - 3?tt- 2 J' m - m 3 K m ) 

+w 3 M m + W 4 (M' m + mN m ) + w 5 

(6.25) 

w 0 K m + Wi(K' m - mJ m ) + w 2 {JZ + - 3 m 2 J' m - m 3 K m ) 


+w 3 N m + W 4 {N' m - mM m ). 

(6.26) 

We solve the equations for M = M m , N = N m , J = J m and K = K m in Auto-07p, 
coupled to a system to find the travelling wave itself, seeking both travelling wave and 
perturbation periodic in Q with period L. Thus we obtain the solution 

h = J m (C) cosmx + K m (0 sinmx = J m (x — Ut) cosmx + K m (x — Ut) sin mx. (6.27) 

Regardless of the relation between the travelling wave period L and the blowing/suction 
wavenumber to, we see that if J m (C) an d K-m(C) are periodic with period L , then the 
function h(x, t) is temporally periodic with period T = U/L , which is the same temporal 
period as the base travelling wave. However, the solution for h.(x, t) is spatially periodic 
only if L/Lp is rational, with period given by the least common multiple of L and Lp. 
A typical set of solutions for H{ £), J m (0 an d K m (C), and also the reconstructed field h, 
are shown in figure m For the travelling wave solution shown in figure EH we find that 
the functions J m (C) aR d K m {Q display little relative variation with £. This means that 
the field h(x,t) is essentially steady. 

Figure [12] shows a comparison between fully nonlinear time-dependent calculations and 
the perturbed travelling wave solution at a value of R small enough that the uniform state 
is unstable to only a single perturbation, and the travelling wave is stable. We obtain 
good agreement between the time-dependent calculations and the asymptotic predictions 


based on (6.10) and (6.27). Both sets of results show a smooth transition as A increases 


from travelling waves with wavelength 60 at A = 0 to almost-steady states at A = 0.18 
with wavelength 20, which is wavelength of the blowing and suction. 

The composite solution shown in figure [Tl] has Lp = 20 and L = 30; the resulting 
perturbed field has spatial period 60. We have only considered corrections up to O(A), 
and in this expansion, the imposed blowing and suction cannot affect the periodicity of 
the underlying travelling wave. However, in fully nonlinear time-dependent simulation, 
even if we begin with initial conditions that are perfectly periodic with period L = 30, 
but force at a different wavelength, such as Lp = 20, nonlinear effects will lead to a 
perturbation at wavelength 60 that can cause the underlying travelling wave to double in 
spatial period, and likely change shape and speed. Eventually we would expect to reach 
the state where the dominant wave has period 60, and the linear perturbation field h is 
a solution to equations where the coefficients Wi are those for the travelling wave with 
wavelength 60. Figure [13] shows time-dependent simulations for forcing at wavelength 20 
with A = 0.02 in a domain of length 60, at a value of R large enough that travelling 
waves with wavelengths 60, 30 and 20 are unstable. When the initial conditions involve 
only modes of wavelength 60, a periodic initial condition is quickly reached. For initial 
conditions of wavelength 30, these modes compete with the wavelength 20 forcing, but 
eventually a wavelength 60 state is indeed achieved. For initial conditions with wavelength 
equal to the initial forcing, we rapidly reach an initial periodic state, with three equal 
waves in the domain. However, after a very long time, noise in the system causes a switch 
to a single wave with wavelength 60, thus yielding the same state regardless of initial 
conditions. 
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Figure 11: Travelling wave H(Q and 0(A) perturbation h(x,t) for 9 = 7r/4, C = 0.05, 
R = 3, calculated with the weighted-residual equations. The initial travelling wave has 
wavelength 30, but blowing and suction is applied with wavelength 20; we show solutions 
in a domain of length 60. The black contours indicate h = 1, and the same colour map is 
used in the first three figures. Nonlinear time-dependent results for A = 0.02 are shown 
in figure 


13 


for a selection of initial conditions. 


7. Initial value problems 

Even in the absence of blowing and suction, a thin liquid film falling along an in¬ 
clined plane can display rich behaviour, including pattern formation, transition to chaos, 
travelling waves, and pulse-like structures. Many of the instabilities that dominate the 
dynamics occur at long wavelength, and so the system is frequently studied using long¬ 
wave models. However, a well-known feature of the Benney model is that the interface 


can exhibit finite-time flow up at Reynolds numbers larger than critical (Pumir et al. 
1983 |Ruyer-Quil fe Mannevillc 2000), which is not replicated in full Navier-Stokes sim¬ 


ulations. The weighted-residual equations were developed partly to avoid this blow-up 
behaviour, and are indeed better behaved that the Benney models. The two models agree 
when the system parameters are close to neutral stability. 

The introduction of forcing in the form of suction boundary conditions introduces 
another potential mechanism for finite-time blow-up, whereby the film thins to zero 
in finite time, as illustrated in figure [4] This thinning occurs even at zero Reynolds 
number, and so arises in both Benney and weighted-residual models. The blow-up need 
not occur at the same wavelength as the forcing function; for example it can be triggered 
by subharmonic-perturbations to a periodic steady state. 











































Figure 12: Contour plots of h(x,t) from nonlinear time-dependent calculations (top row) 
and using the small-A asymptotic solutions for perturbed travelling waves (bottom row) 
for R = 1.7, C = 0.05, 6 = 7t/4, in a domain of length 60 with suction wavelength 
Lp = 20. The colour map is scaled to the maximum and minimum value in each column. 
Here the travelling wave is stable in the absence of suction. Time-dependent results for 
a larger R, for which the uniform film solution is unstable to multiple perturbations, are 


shown in figure 13 


Blow-up is generally avoided if R is not too large, and blowing and suction is applied 
with sufficiently small amplitude, but the film dynamics can still interact with the im¬ 
posed suction. In the absence of suction, the system exhibits periodic behaviour in the 
form of travelling waves, which propagate at a constant speed without changing form. If 
non-uniform suction is imposed via a non-constant function F{x ), which remains fixed 
with respect to the frame of the wall, waves must change in form as they propagate, 
and so time-periodic behaviour manifests as nonlinear limit cycles which depend on both 
x and t. Figure [12] shows initial value calculations in which the unforced system ex¬ 
hibits stable travelling waves. The imposed suction need not have the same wavelength 
as the travelling wave, and in figure three periods of suction fit within the discre- 
tised domain. The initial value calculations show a smooth transition as A increases from 
travelling waves at A = 0 to an essentially steady state at A = 0.18. At small A, the 
blowing and suction causes slight perturbations to the travelling wave, while at larger 
A, small-amplitude disturbances propagate over a large-amplitude non-uniform steady 
state. The underlying flow field retains its dominant down-slope direction, and so there 
is still a non-zero perturbation wave speed even when the interface shape is steady. 

The applied suction can also lead to states which are neither steady nor limit cycles. 
Figure [14] shows the result of initial value calculations conducted at parameters near to 
the stable ‘island’ shown in figure [9] for L = 10. This is a region of parameter space with 
steady solutions that are stable to perturbations of all wavelengths, but is unusual in 
that solutions lose stability if either A or R is decreased, so that sufficiently large inertia 
is required to maintain stability. In case shown in figure [l4]( a), with A = 0.25 and R = 4, 
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Figure 13: Time-dependent simulations in a domain of length 60, with F = 
0.02cos(27ra:/20), with R = 3, C = 0.05 and 9 = 7r/4, and varying initial conditions: we 
set h(x, 0) = 1 + 0.1 cos(27rna;/60) and q(x, 0) = 2/3 + 0.2cos(27rn:r/60), with n = 1, 2,3 
in (a-c) respectively. These initial conditions correspond to neither travelling waves nor 
steady states, so nonlinear evolution is required if the system is to reach a periodic state. 
These plots shows the single contour h(x,t) = 1, with h > 1 in the shaded region. The 
same periodic state is reached eventually, regardless of the initial conditions. 


the height field displays aperiodic behaviour. At each instant, six peaks corresponding to 
the six periods of the suction function are visible, but waves propagate over these peaks 
without a clear structure. We can increase the suction amplitude to A = 0.45 to obtain 
steady, stable solution, shown in figure [l4| 6). Decreasing the Reynolds number to R = 2 
means that the steady state again loses stability, but in this case to a wavelength 60 
travelling wave which propagates over the steady solution in a regular manner, as shown 
in figure [life). 


8. Conclusion 

In this paper, we considered thin-film flow down an inclined plane modified so that fluid 
is injected and withdrawn through the wall according to an arbitrary function F[x , t). We 
derived and studied two long-wave models, based on the first-order Benney and weighted- 
residual formulations, for this system. If the average layer height is conserved, then F 
must have zero mean in space. We then specialised to the case where A is a single steady 
Fourier mode, F = Acosmx , and investigated the form, bifurcations and linear stability 
of steady states, as well as a range of fully-nonlinear time-dependent behaviours. 
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(c) A = 0.45, R = 2 



Figure 14: Nonlinear time-dependent calculations for L = 10, C = 0.05, 9 = 7r/4. The 
simulations are conducted in a domain of length 6 L, and the initial conditions for h 
and q includes a small perturbation proportional to sin(27ra;/(6L)). The stable solution 
shown in ( b ) can be accessed by either increasing the Reynolds number from (c), or the 
amplitude of blowing and suction from (a). 


Any steady states subject to non-uniform blowing and suction F must themselves be 
non-uniform. We calculated the interface shape for small A, and showed that this differs 
between the two models when R ^ 0, but that the results agree in the long-wave limit. 
There is a phase difference between the shape of the blowing and suction and interface 
height, which originates from the preferred down-slope direction of the base flow. 

As A increases, the interface becomes increasingly non-uniform, and the nonlinear 
terms in the governing equations become increasingly important. Two notable transitions 
occur as A increases: firstly all steady states feature regions with negative volume-flux 
q, and secondly steady solutions disappear altogether for sufficiently large A. Time- 
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dependent evolution beyond the existence of steady states shows that the film can dry in 
finite time, as the fixed rate of fluid removal via F is not matched by the supply of fluid. 

We analysed the behaviour of the streamfunction in order to interpret steady solutions 
with negative q. We found that imposing non-zero, non-uniform F leads to a series of 
isolated stagnation points along the wall, wherever F(x) = 0. These stagnation points 
are either saddle points or local extrema in the streamfunction, depending on the sign 
of F'(x) and q(x,t) at these points. Connections between the stagnation points allow 
us to demarcate the boundary between fluid which enters and leaves through the wall, 
and ‘propagating’ fluid which never reaches the wall. At A = 0, all fluid propagates 
down the wall, but the height of the recirculating layer near the wall increases with A. 
For large enough A, all steady solutions have negative q for some x. For those x with 
q{x) < 0, all flow is directed up the slope, against gravity; this flow reversal means that 
the propagating layer disappears. 

Steady solutions need not have the same spatial periods as the imposed suction, and 
subharmonic steady solutions can emerge via pitchfork bifurcations from the harmonic 
solution branches. In our calculations, we have only found subcritical pitchfork bifur¬ 
cations, leading to unstable subharmonic steady states. Unstable subharmonic steady 


states also occur in thin-film flow with periodic topography (Tseluiko et al. 2013), but 


it is possible that subharmonic steady states may be stable for other parameter values. 
Initial-value calculations starting from near the periodic steady states sometimes show 
subharmonic drying, where the film dries at any one of the local minima of h. 

In the absence of suction, the primary instability of the uniform state is to eigenfunc¬ 
tions with complex eigenvalues, leading to Hopf bifurcations which do not appear in the 
bifurcation structure for steady solutions. If F is spatially periodic, then the steady states 
are also periodic, as are the coefficients of the linear equation governing the evolution of 
small perturbations. This means that perturbations must propagate through a heteroge¬ 
neous environment, and also that the eigenmodes are no longer simply exponentials of 
the form exp (ikx + Xt). Instead, according to the Floquet-Blocli form, the eigenmodes 
can be written as g(x) exp (ikx + At), where g(x) is a complex-valued function, spatially 
periodic with the same period as the base state. We calculated linear stability numeri¬ 
cally for two classes of perturbations: firstly those that are periodic with the same period 
as the blowing and suction, and secondly perturbations of any wavelength. We found 
that introducing spatially-periodic suction with amplitude A can increase or decrease 
the critical Reynolds number for linear stability to perturbations of both classes. Due to 
the symmetries of the system, the correction to the eigenvalue, and hence to the critical 
Reynolds number, is even in A. 

For small A, we calculated the critical Reynolds number to 0(A 2 ) for perturbations of 
arbitrary wavenumber. Imposing blowing and suction at very long wavelength increases 
the critical Reynolds number for both classes of perturbations when 9 < 7t/2. Due to the 
quadratic dependence of the critical R on perturbation wavenumber, the first eigenmodes 
to become unstable as R increases are those with infinite wavelength, both when A = 0 
and for small A, In contrast, for the case of flow underneath an inclined plane, where 
9 > 7t/2, small-amplitude long-wave blowing and suction always reduces the critical 
Reynolds number, and so has a destabilising effect on the flow. However, as the small 
A expansion is valid only close to the critical Reynolds number, which is negative for 
9 > 7t/2, the predicted dependence of stability properties on A may not hold for physically 
realisable flows. 

Determination of the largest R that can be stabilised in an infinite domain by imposing 
steady suction is difficult as it is not clear that long-wavelength waves are always the most 
unstable, and both the steady solutions and their corresponding linear stability operators 
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must be computed at finite A. We find that the weighted-residual model predicts large 
increases in the critical Reynolds number for stability to perturbations of all wavelengths, 
from 1.25 to 4 or more, if blowing and suction is introduced with wavelengths of around 
200. However, the imposed suction appears to become less effective at stabilising the flow 
if applied with very long wavelength. 

In the absence of suction, the system can support travelling waves which propagate at 
a constant speed U without changing form. Periodic travelling waves are periodic with 
respect to both space and time, but their dependence on x and t can be expressed in 
terms of a single variable £ = x — Ut. Imposing blowing and suction as F(x) introduces 
heterogeneity to the system, and so we expect travelling waves to become limit cycles that 
are periodic in both x and t, with explicit dependence on both variables. We calculated the 
O(A) corrections to the travelling wave, and showed that the equations are closely related 
to those governing linear stability of the travelling wave. We demonstrated that limit 
cycles at O(A) can be decomposed by rewriting the equations in terms of the independent 
variables x and £. The resulting solution for h(x,t) and q(x,t) is periodic with respect to 
time with the same time period as the underlying travelling wave, but is spatially periodic 
only if the ratio between the wavelength of the blowing and suction and the wavelength of 
the travelling wave is rational. We showed that the predictions of the O(A) perturbation 
analysis are in good agreement with the results of fully nonlinear simulations. We also 
conducted initial-value calculations to investigate competition between the period of the 
travelling wave and of the imposed suction. We found that the same state is eventually 
reached, but this competition can persist over a large number of cycles. 

We observed in § |2.3| that the derivation of our equations is unchanged if the blowing 
and suction profile additionally depends on time with sufficiently long timescale. Allowing 
time-dependence of the blowing and suction leads to the possibility of using F to explicitly 
control the flow, and thus to stabilise otherwise unstable states by acting in response 
to the development of perturbations. For example, with perfect implicit knowledge of 
the system, we could choose F(x,t) in order to drive the system towards a particular 
steady state. However, the additional independent degree of freedom that emerges in 
the weighted-residual model may be important when exploring control of the system, 
as we cannot choose F(x,t) to simultaneously specify h(x,t ) and q{x,t). Under the 
assumption of small deviations from a uniform state, with small F, both the weighted- 
residual equations and the Benney equations reduce via a weakly nonlinear analysis to 
a forced version of the Kuramoto-Sivanshinsky equations, for which optimal controls 
have been successfully applied ( Gomes et aL|2015 1. Work is in progress (Thompson et al. 
2015) to explore the effectiveness of feedback control strategies based on the nonlinear 


long-wave models derived in this paper. 
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Appendix A. Derivation of long-wave equations 

We first rescale the governing equations (2.1) to (2.7) according to 

X = 6x, T = St, v = 6w, C = 8 2 C , (Al) 

where we are interested in the long-wave limit { « 1. We will consider two different 
scalings for F: firstly F = Sf, and also F = S 2 f. If F is unsteady, we require for 
consistency that fj~ = 0(1) or smaller. The full set of equations becomes 

RS (ut + uux + wu y ) = -5px + 2 + 5 2 uxx + u vy , 

R5 2 (wt + wwx + wvjy) = — p y — 2 cot 8 + 8 3 wxx + 8w y 
Ux + Wy = 0 . 

The boundary conditions at y = 0 are 

u = 0, w = f(X,T), 

and at y = h we have 

(S 2 w x + Uy) (l - 5 2 h 2 x ) + 2 5h x ( Sw v - Su x ) = 0, 

-- [Swy + 5 3 uxh\ - Sh x (6 2 w\ + u y )) = ——- ^ xx — (A 7) 

1 + 5 2 h 2 x y v x K vn 0(1 + S 2 h 2 x ) 3 / 2 


u yy> 


(A 2) 
(A3) 
(A 4) 

(A 5) 

(A 6) 


P~Pa- 


The rescaled kinematic equation (2.8) is 


hr ~ f(X,T) +q x = 0, q=[ udy. 

Jo 

As / is known, we need only an expression for q to close the system. 


(AS) 


A.l. Benney equation 

To obtain the Benney equation for q as a function of h, we expand u, w, p and q in 
powers of 5, while assuming that h is an 0(1) quantity: 

u = uq + Sui + 0{8 2 ), w = wq + dioi + 0(8 2 ) (A 9) 

p = Po + Spi + 0(S 2 ), q = qo + 5q± + 0(S 2 ). 


(A 10) 


The leading-order solution of (A 2) to (A 7) for small 6 is 


h.xx 


u 0 = y(2h-y), w 0 = f{X,T) - yhx, Po = Pa -— + 2(h - y) cot 8. (All) 


C 


The flux at this order is 


f h , 2 h 3 

qo= u 0 dy= —, 


(A 12) 
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leading to the evolution equation 

h T - f(X,T) + 2h 2 h x = 0(6). 


(A 13) 


We must calculate iq in order to obtain the 0(6 ) correction to q. After integrating the 
0(5) part of (A2) twice, and applying boundary conditions at this order, we find 


ui = y -2 h) + R 


where we have defined 


(hT-f)(*j-h 2 y) + 2hhx (V ' 


- h 3 y + hy(y - 2 h)f 


p = 2h cot 0 — 


ixx 

C 


The first-order correction to the flux can now be calculated as 

3 h 6 hx 
10 


Qi = mi d y = - 


h 3 p 


'x 


R - 


5 hrh 4 
12 




(A 14) 
(A 15) 

(A 16) 


We eliminate hx from this expression by using (A 13), and thus obtain the Benney 
equation for q: 


q = q 0 + 6qi + 0(6 2 ) = ^-+ RS 


( 8 h e h x 2 h A f(X,T) 

V 15 3 


+ 0(8 2 ). (A 17) 


Alternatively, if / == 6f , then / still appears via the mass conservation (A8), but is 
absorbed into the 0(5 2 ) error term in (A 17). 

A.2. First-order weighted-residual equations 

The derivation of the first-order weighted-residual equations closely follows the original 
derivation presented by Ruyer-Quil & Manneville (2000). Imposing suction at the wall 
does not affect the boundary conditions on u, and so the basis functions for u described 
by |Ruyer-Quil fe Manneville] yield the first-order equation without difficulty. 

In contrast to the derivation of the Benney equations, here we use 6 as an ordering 
parameter, rather than directly expanding variables with respect to <5. For the first-order 
weighted-residual equations, we retain terms up to and including 0(5) in the equations, 
and so the momentum equations yield 

R (Sut + Suux + Swu.y) = —Spx + 2 + u yy + 0(5 2 ) 


and 


0 = — Py — 2 COt 9 + 5w yy + O (f^), 


(A 18) 


(A 19) 


while the mass conservation equation (A 4) and the boundary conditions on the wall 
(A 5) are unchanged. The normal and tangential components of the dynamic boundary 

hxx 


condition become 

0 = u y + 0(5 2 ) 


p = p a + 25ux -A—b 0(5 2 ) at y = h. 

C 


(A 20) 


(A 21) 


We can integrate (A 19) with respect to y. and apply (A 20) to obtain 

p = p a + 2 (h — y) cot 6 — ~ $ u x — 2 5hxu y . 

We then substitute (A21) into ( A 18| ) and discard terms smaller than 0(5 ), leaving 

R6 (u T + uu x + wu v ) = -5p x + 2 + u yy , (A 22) 
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which is coupled to the mass conservation equation (A4), and subject to the boundary 
conditions u = 0 and w = f(X,T) at y = 0, and u y = 0 at y = h(X,T). 

Following Ruyer-Quil & Manneville, we posit an expansion for u in terms of basis 
functions <f>j satisfying no-slip on the wall and zero tangential stress on the interface: 


= '}2a j {X, T )cf> j (y), (f>j(z) = z J+1 - 


j + 1 

3 + 2 


y0+ 2 


y = 


h{X, T )' 


(A 23) 


If 6 = 0, the only non-zero a n is a o, and for small S, we find that a o = 0(1), while 
a n = 0(5) or smaller for n > 1. We can calculate the coefficients aj by using the tpj as 
test functions in (A22), but can obtain the leading order equations by use of cf>o only. 
Correct to 0(5) we can write the weak form of (A22) as 

nh nh nh 

R5 (fn(y) (u 0 T + u 0 u 0 x + W 0 u 0y ) dy = (~5px + 2) 4>n(y)dy + <t> n (y)u yy dy, 

J o Jo Jo 

(A 24) 

where 

u 0 (X,y,T)= 3 ^My), w 0 (X,y,T) = f(X,T)-J\ox(X, y ',T)dy'. (A 25) 

After setting n = 0 and repeated integration by parts, (A24) yields 

18 q 2 hx 


( 2 23 qh T 

M| 5 ,T- 40T 


111 qqx_ _ 

280 h 35 h 2 


+ ¥)=(-^+2 (A 26) 


8 h 


We can eliminate hx from (A26) by using (A8), and thus obtain 

hqf(X,T) 


, ^ D^J .2 2h3 fih 3 px 

q+-R5h q T = — - — 

5 3 3 


nf , 18 2l 34 
+ R5 ( — q h x - + 


+ 0(5 2 ). (A27) 


The two equations (A 8) and (A 27) form a closed system for h and q without requiring 
calculation of a n for n > 1. However, the higher coefficients are required to fully determine 
the velocity field within the fluid layer. 

As in the case of the Benney equations, if we take / = 5f, the mass conservation equa¬ 
tion ( |A8[ ) is unchanged, but the term involving / in (A27) disappears as it is absorbed 
into the 0(5 2 ) error term. 




















